
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON /ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON/HLP/NBLDT
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /MESHSP/SPHALF,SKR,SKT
      COMMON /XRTSAV/XSAVE(50,100,3),RSAVE(50,100,3)
      COMMON/TTSAVE/TTSAV(50,50,3),TRSAV(50,50),TXSAV(50,50),
     *   TPSAV(50,100),TSSAV(50,100)
      COMMON/FILEWR/IWRITE
C
      READ(5,500)(TITLE(I),I=1,80)
500   FORMAT(80A1)
      READ(5,510)NBLADE
510   FORMAT(5I5)
      NBLDT=NBLADE
C
C
      READ(5,510)NX,NTH,NR
      READ(5,510)N1,N2,KUTTA,IEXPLE,IEXPTE
      IF(NX.LE.100)GO TO 15
      WRITE(6,325)NX
      GO TO 1000
15    IF(NTH.LE.20.OR.NR.LE.20)GO TO 25
      IF(NTH.GT.50.OR.NR.GT.50)GO TO 20
      WRITE(6,326)NTH,NR
      GO TO 25
20    WRITE(6,327)NTH,NR
      GO TO 1000
25    CONTINUE
      IF(N2.LE.N1)GO TO 35
      IF(N1.LE.0)GO TO 35
      IF(N2.GE.NX)GO TO 35
      GO TO 40
35    WRITE(6,328)NTH,NR
      GO TO 1000
40    IF(IEXPLE.EQ.0) GO TO 50
      IF(IEXPLE.LE.N1 .OR. IEXPLE.GE.N2) GO TO 45
      IF(IEXPTE.LE.N2.OR.IEXPTE.LE.IEXPLE)GO TO 45
      GO TO 50
45    WRITE(6,329)IEXPLE,IEXPTE
      GO TO 1000
50    CONTINUE
      DX= 2.0/(N2-N1)
      DR=1.0/(NR-2)
      DTH=1.0/(NTH-2)
      READ(5,520)SPHALF
520   FORMAT(3F10.5)
      SKR=0.0
      SKT=0.0
      IF(SPHALF.NE.0.)SPHALF=0.5
C
      IF(SPHALF.EQ.0.)DR=1.0/(NR-1)
      IF(SPHALF.EQ.0.)DTH=1./(NTH-1)
      READ(5,510)IAN
      IF(IAN.NE.1)IAN=0
C
C
      WRITE(6,315)(TITLE(I),I=1,80),NBLADE
315   FORMAT(' 3D AXIAL COMPRESSOR GRID CODE.',/,
     11X,80A1,/,' THE NUMBER OF BLADES IS ',I5)
      WRITE(6,324)SPHALF
324   FORMAT(' HALF MESH SPACING CONSTANT IS ',F5.3)
      WRITE(6,318)NX,NTH,NR
318   FORMAT(' GRID POINTS IN X DIRECTION = ',I5/
     *' IN THETA DIRECTION = ',I5,' IN R DIRECTION =',I5)
      WRITE(6,319)N1,N2,KUTTA
319   FORMAT(' THE BLADE LEADING EDGE AND TRAILING EDGE ARE'
     *,' AT AXIAL GRID LINES ',2I5/' KUTTA POINT IS',
     *' AT GRID LINE ',I5)
      WRITE(6,322)IEXPLE,IEXPTE
322   FORMAT(' THE FIRST POINT ON THE DAMPER IS ',I5/
     *' THE LAST POINT ON THE DAMPER IS ',I5)
      IF(IAN.EQ.1)WRITE(6,334)
      IF(IAN.NE.1)WRITE(6,337)
C
C
      IER=0
      READ(5,510)IWRITE,IDMP
      WRITE(6,335)IWRITE,IDMP
C
      WRITE(IWRITE,315)(TITLE(I),I=1,80),NBLADE
      WRITE(IWRITE,318)NX,NTH,NR
      WRITE(IWRITE,319)N1,N2,KUTTA
      WRITE(IWRITE,322)IEXPLE,IEXPTE
C
      CALL XRGRID(IWRITE,IDMP,IER,IAN)
      WRITE(6,330)IER
      IF(IER.EQ.1)GO TO 1000
      CALL CALT(IWRITE,IDMP,IER)
      WRITE(6,331)IER
      CALL BLDT(IWRITE,IDMP,IAN)
      WRITE(6,332)
      CALL BUILD(IWRITE)
      WRITE(6,333)
      CALL SVGRID
      WRITE(6,336)
C
C
1000  CONTINUE
325   FORMAT(' NX=',I5,' IS GREATER THAN NXMAX.')
326   FORMAT(' NTH OR NR',2I5,' IS TOO LARGE FOR INVISCID CODE',
     *         ' BUT WITHIN RANGE FOR METRIC CODE')
327   FORMAT(' NTH OR NR',2I5,'IS TOO LARGE FOR METRIC CODE')
328   FORMAT(' N1 OR N2(LEADING OR TRAILING EDGE POSITION)',
     *' SPECIFIED INCORRECTLY.',2I5)
329   FORMAT(' IEXPLE OR IEXPTE',2I5,' IS INCORRECT')
330   FORMAT(' SUBROUTINE GRID FINISHED.IER=',I5)
331   FORMAT(' SUBROUTINE CALT FINISHED.IER=',I5)
332   FORMAT(' SUBROUTINE BLDT FINISHED')
333   FORMAT(' SUBROUTINE BUILD FINISHED')
334   FORMAT(' ANALYTICAL METRICS USED')
337   FORMAT(' NUMERICAL METRICS USED')
335   FORMAT(' OUTPUT FILE NUMBER IS',I5/
     *' DEBUG DUMP PARAMETER IS ',I5)
336   FORMAT(' SUBROUTINE SVGRID FINISHED')
C
C
C
      STOP
      END
      SUBROUTINE ARC(N,X,R,T,AL)
      REAL X(N),R(N),T(N),AL(N)
      AL(1)=0.
      DO 100 J=2,N
      RS= R(J)
      J1=J-1
      R2=RS**2
      S1=(RS-R(J1))**2
      S1=(X(J)-X(J1))**2 +S1
      AL(J)=SQRT(S1)+AL(J1)
100   CONTINUE
      RETURN
      END
      SUBROUTINE BLDT(IWRITE,IDMP,IAN)
C
      REAL DTPDR(50),DTSDR(50),RES(2,50),RTAU(2,50)
      COMMON /MESHSP/SPHALF,SKR,SKT
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/HLP/NBLDT
      REAL TR(50,50),TX(50,50),TT(50,50),TR1(50),TH(50)
      REAL CUR(50,2,2),COSN(50,3,2)
      COMMON /ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON /XRTSAV/XSAVE(50,100,3),RSAVE(50,100,3)
      COMMON/TTSAVE/TTSAV(50,50,3),TRSAV(50,50),TXSAV(50,50),
     *   TPSAV(50,100),TSSAV(50,100)
      REAL CR(2,50),CT(2,50),CZ(2,50)
      REAL X(50),R(50),TS(50),TP(50),DTPDZ(50),DTSDZ(50)
      NR1=NR-1
      NTHA=NTH-2
      IWR=IWRITE
      NR2=NR-2
      NTH1=NTH-1
      N1M1 = N1-1
      N2P1 = N2+1
      PI=3.14159265
      PPSAV=2.*PI/NBLADE
      PPSAV=1./PPSAV
      REWIND 1
      READ(9)X,R,TS,TP,DTPDZ,DTSDZ,DTPDR,DTSDR
      REWIND 9
      CDZA=0.
      DO 1004 L=1,NR
      CDZA=CDZA+(DTPDZ(L)+DTSDZ(L))*0.5
1004  CONTINUE
      CDZA=CDZA/NR
      CDRA=0.
      DO 4 L=1,NR
      DTHICK=(2.*PI/FLOAT(NBLDT)-(TS(L)-TP(L)))/2.
C     WRITE(6,*)TS(L),TP(L)
      CDZS=(DTPDZ(L)+DTSDZ(L))*0.5
      CDRS=(DTPDR(L)+DTSDR(L))*0.5
      CPL=TP(L) - DTHICK
      CSL=TS(L) + DTHICK
      DO 4 J=1,N1M1
      TTMP=(N1-1)/J
      TTMP=EXP(-TTMP)
      CDZ=CDZA+(CDZS-CDZA)*TTMP
      CDR=CDRA+(CDRS-CDRA)*TTMP
      CZZ=XSAVE(L,N1,1)-XSAVE(L,J,1)
      CRR=RSAVE(L,J,1)-RSAVE(L,N1,1)
      TPSAV(L,J)=CPL-CDZ*CZZ+CDR*CRR
      TSSAV(L,J)=CSL-CDZ*CZZ+CDR*CRR
4     CONTINUE
      DO 5 J=N1,N2
      READ(4)TS,TP
      DO 5 L=1,NR
C     WRITE(6,*)TS(L),TP(L)
      TPSAV(L,J)=TP(L)
5     TSSAV(L,J)=TS(L)
      REWIND 4
      DO 6 J=N1,N2
      READ(9)X,R,TS,TP,DTPDZ,DTSDZ,DTPDR,DTSDR
      READ(9)CR,CT,CZ
6     READ(9)RES,RTAU
      REWIND 9
      CDZA=0.
      CDRA=0.
      DO 1007 L=1,NR
      CDZA=CDZA+(DTPDZ(L)+DTSDZ(L))*0.5
1007  CONTINUE
      CDZA=CDZA/NR
      DO 7 L=1,NR
      DTHICK=(2.*PI/FLOAT(NBLDT)-(TS(L)-TP(L)))/2.
      CDZS=(DTPDZ(L)+DTSDZ(L))*0.5
      CDRS=(DTPDR(L)+DTSDR(L))*0.5
      CPL=TP(L) - DTHICK
      CSL=TS(L) + DTHICK
      DO 7 J=N2P1,NX
      TTMP=(NX-N2+1)/(NX-J+1)
      TTMP=EXP(-TTMP)
      CDZ=CDZA+(CDZS-CDZA)*TTMP
      CDR=CDRA+(CDRS-CDRA)*TTMP
      CZZ=XSAVE(L,J,1)-XSAVE(L,N2,1)
      CRR=RSAVE(L,J,1)-RSAVE(L,N2,1)
      TPSAV(L,J)=CPL+CDZ*CZZ+CDR*CRR
      TSSAV(L,J)=CSL+CDZ*CZZ+CDR*CRR
7     CONTINUE
      DO 20 J=1,NX
      DO 18 L=1,NR
      DO 17 K=1,NTH
      TR1(K)=K
17    TH(K)=TPSAV(L,J)+(TSSAV(L,J)-TPSAV(L,J))*(K-1)/(NTH-1)
      IF(SPHALF.EQ.0.) GOTO 32
      DELTHE=(TSSAV(L,J)-TPSAV(L,J))/(NTH-2)
      TH(1)=TPSAV(L,J)-0.5*DELTHE
      DO 31 K=2,NTH
      TH(K)=TPSAV(L,J)+DELTHE*(FLOAT(K)-1.5)
31    CONTINUE
32    CONTINUE
      IF(SKT.NE.0.)CALL RECLUS(SKT,NTH,TR1,TH)
      DO 19 K=1,NTH
19    TTSAV(K,L,1)=TH(K)
18    CONTINUE
      WRITE(1)((TTSAV(K,L,1),K=1,NTH),L=1,NR)
20    CONTINUE
      REWIND 1
      READ(1)((TTSAV(K,L,2),K=1,NTH),L=1,NR)
      READ(1)((TTSAV(K,L,3),K=1,NTH),L=1,NR)
      DO 100 J=1,NX
      IF(J.EQ.2)GOTO 23
      IF(J.EQ.NX)GOTO 23
      DO 22 K=1,NTH
      DO 22 L=1,NR
      TTSAV(K,L,1)=TTSAV(K,L,2)
22    TTSAV(K,L,2)=TTSAV(K,L,3)
      READ(1)((TTSAV(K,L,3),K=1,NTH),L=1,NR)
23    CALL YYM(J)
      DO 26 L=1,NR
      DO 24 K=1,2
      CUR(L,K,1)=0.
24    CUR(L,K,2)=0.
      DO 25 K=1,3
      COSN(L,K,1)=0.
25    COSN(L,K,2)=0.
26    CONTINUE
      IF(J.EQ.1)READ(9)X,R,TS,TP,DTPDZ,DTSDZ,DTPDR,DTSDR
      IF(J.EQ.1)REWIND 9
      IF(J.LT.N1) GO TO 28
      IF(J.GT.N2)GO TO 28
      READ(9)X,R,TS,TP,DTPDZ,DTSDZ,DTPDR,DTSDR
      READ(9)CR,CT,CZ
      READ(9)RES,RTAU
28    CONTINUE
      JP=J-N1+1
      DO 40 L=1,NR
      IF(IAN.EQ.0)GO TO 43
      IF(J.GE.N1 .AND. J.LE.N2)GO TO 41
      TP(L)=0.
      TS(L)=1./PPSAV
      T1=(DTSDR(L)+DTPDR(L))*0.5
      DTSDR(L)=T1
      DTPDR(L)=T1
      T1=(DTSDZ(L)+DTPDZ(L))*0.5
      DTPDZ(L)=T1
      DTSDZ(L)=T1
41    CONTINUE
      T1=TS(L)-TP(L)
      T2=1./T1**2
      OTH=T1/(NTH-2)
      T3=DTSDR(L)-DTPDR(L)
      T4=DTPDR(L)*T1
      T5=DTSDZ(L)-DTPDZ(L)
      T6=DTPDZ(L)*T1
      L1=L
      THETA=TP(L)-OTH/2.
      DO 42 K=1,NTH
      P=THETA-TP(L)
      P1=P*T3+T4
      TRSAV(K,L)=-P1*T2
      P1=P*T5 +T6
      TXSAV(K,L)=-P1*T2
      TPSAV(K,L)=1./T1
      THETA=THETA+OTH
42    CONTINUE
43    CONTINUE
      L1=L
      DO 30 K=1,NTH
      TT(K,L)=TPSAV(K,L)
      TR(K,L)=TRSAV(K,L)
      TX(K,L)=TXSAV(K,L)
30    CONTINUE
      IF(J.LT.N1)GO TO 40
      IF(J.GT.N2)GO TO 40
C       FIND GRAD CAPTHETA COMPONENTS AT PRESSURE AND SUCTION SURFACES
C
C       PRESSURE SURFACE
C
      COSN(L1,1,1)=CR(1,L)
      COSN(L1,2,1)=CT(1,L)
      COSN(L1,3,1)=CZ(1,L)
      CUR(L1,1,1)=RES(1,L)
      CUR(L1,2,1)=RTAU(1,L)
C
C       SUCTION SURFACE
C
      COSN(L1,1,2)=CR(2,L)
      COSN(L1,2,2)=CT(2,L)
      COSN(L1,3,2)=CZ(2,L)
      CUR(L1,1,2)=RES(2,L)
      CUR(L1,2,2)=RTAU(2,L)
40    CONTINUE
      WRITE(IWR,700)J
      WRITE(IWR,701)
      DO 710 KX=1,NTH
      DO 710 LX=1,NR
      IF(KX.NE.1.AND.KX.NE.NTH.AND.MOD(KX,5).NE.0)GO TO 710
      IF(LX.NE.1.AND.LX.NE.NR .AND.MOD(LX,5).NE.0)GO TO 710
      WRITE(IWR,702)KX,LX,TT(KX,LX),TR(KX,LX),TX(KX,LX)
710   CONTINUE
700   FORMAT(///10X,'THET.DER. FOR J=',I4)
701   FORMAT(/6X,'K',3X,'L',5X,'TT ',7X,'TR ',7X,'TX ')
702   FORMAT(5X,I3,1X,I3,1X,E9.3,1X,E9.3,1X,E9.3)
      IF(J.LT.N1)GO TO 50
      IF(J.GT.N2)GO TO 50
      WRITE(IWR,600)J
600   FORMAT(//,'  NORMAL VECTORS AND RADII OF CURVATURE FOR AXIAL STATI
     1ON NO. ',I5)
      WRITE(IWR,601)
601   FORMAT(1X,T4,3('L  COSNR  COSNT  COSNZ    RES   RTAU',6X))
      WRITE(IWR,604)(L,(COSN(L,KX,1),KX=1,3),RES(1,L),
     *   RTAU(1,L),L=1,NR)
      WRITE(IWR,605)(L,(COSN(L,KX,2),KX=1,3),RES(2,L),
     *   RTAU(2,L),L=1,NR)
604   FORMAT('  PRESSURE SIDE NORMAL VECTOR '
     *  /(3(I4,5F7.3,3X)))
605   FORMAT('  SUCTION SIDE NORMAL VECTOR'
     *  /(3(I4,5F7.3,3X)))
50    CONTINUE
      WRITE(2)((TT(K,L),K=1,NTH),L=1,NR),
     *   ((TR(K,L),TX(K,L),K=1,NTH),L=1,NR),
     *   (((CUR(L,KX,KY),KX=1,2),(COSN(L,KX,KY),KX=1,3),KY=1,2),L=1,NR)
100   CONTINUE
      REWIND 2
      REWIND 9
      RETURN
      END
      SUBROUTINE BUILD(IWRITE)
      REAL R(50),XX(50),XR(50),RX(50),MACH
      COMMON /PROP/GAMMA,W,G1
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/DAMPPS/IDLOW,IDHIGH
      REAL RVD(4)
      REAL JAC(50),DTDT(50,50)
      REAL RV(4),RR(50)
      REAL ZS(50)
      REAL RRS(50),RXS(50),XXS(50),XRS(50)
      IWR=IWRITE
      NR1=NR-1
      DO 50 J=1,NX
      MACH=0.50
      UR=0.
      UTH=0.
      IDAMP=0
      JJ=J-1
      READ(2)((DTDT(K,L),K=1,NTH),L=1,NR)
      NR1=NR-1
      DO 10 L=1,NR
      READ(3)ZS(L),R(L),RRS(L),RXS(L),XXS(L),XRS(L)
10    CONTINUE
      X=(J-1)*DX
      READ(3)DRHDX,D2RHDX
      READ(3)DRTDX,D2RTDX
      CALL SETR(DRHDX,DRTDX,RV,RTHUB,RTTIP,D2RHDX,D2RTDX)
      WRITE(IWR,601)J,DRHDX,RTHUB,DRTDX,RTTIP
601   FORMAT(/,T5,'HUB SLOPE    HUB CURV.    TIP SLOPE    TIP CURV.',
     1' FOR AXIAL STATION NO. ',I5,/,4(3X,F10.5))
      IF(J.LT.IEXPLE)GO TO 18
      IF(J.GT.IEXPTE)GO TO 18
      IDAMP=1
      READ(3)DRDAML,D2DAML,DRDAMU,D2DAMU
      CALL SETR(DRDAMU,DRDAML,RVD,RTUP,RTLOW,D2DAMU,D2DAML)
      DPR=R(IDHIGH)-R(IDLOW)
      RA=R(IDHIGH)
      RB=R(IDLOW)
18    CONTINUE
      DO 20 L=1,NR
      XX(L)=XXS(L)
      RX(L)=RXS(L)
      XR(L)=XRS(L)
      RR(L)=RRS(L)
20    CONTINUE
      IF(R(1).GT.0)GO TO 22
      R(1)=R(2)
      XX(1)=XX(2)
      RX(1)=RX(2)
      XR(1)=XR(2)
      RR(1)=RR(2)
22    CONTINUE
      DO 24 L=1,NR
      D=RR(L)*XX(L)-RX(L)*XR(L)
      LL=L-1
      IF(LL.LT.1)LL=1
      IF(LL.EQ.NR1)LL=NR1-1
      D=D*DTDT(2,L)
      JAC(L)=D
24    CONTINUE
      WRITE(IWR,602) (L,JAC(L),L=1,NR)
602   FORMAT( ' JACOBIANS FOR RADIAL GRID LOCATIONS.'/(9(I3,F10.5)))
      WRITE(4)(R(L),L=1,NR),(XX(L),L=1,NR),(XR(L),L=1,NR),
     *(RX(L),L=1,NR),(RR(L),L=1,NR)
      DMM=IDAMP
      WRITE(4)(RV(I),I=1,4),RTHUB,RTTIP,DMM,RA,RB,(RVD(I),I=1,4),
     *   RTLOW,RTUP
      WRITE(4)(JAC(L),L=1,NR)
50    CONTINUE
      REWIND 3
      RETURN
      END
      SUBROUTINE CALDR(IDMP)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      REAL TS(50),TP(50),TI(50),D2S(50),D2P(50),RS(50),TSS(50)
      REAL X(50),R(50),SS(50),SP(50),DTP(50),DTS(50),XS(50),TPS(50)
      JPMAX=N2-N1+1
      IBUG=0
      IF(IDMP.LT.0  .OR. IDMP.EQ.4)IBUG=1
      DO 100 L=1,NR
      DO 10 J=N1,N2
      READ(2)X,R
      READ(4)TS,TP
      JP=J-N1+1
      XS(JP)=X(L)
      RS(JP)=R(L)
      TPS(JP)=TP(L)
      TSS(JP)=TS(L)
10    CONTINUE
      REWIND 2
      REWIND 4
C
C
C
      CALL ARC(JPMAX,XS,RS,TPS,SP)
      CALL ARC(JPMAX,XS,RS,TSS,SS)
C
C
      CALL INPR(JPMAX,SS,TSS,DTS,D2S,JPMAX,SS,XS)
      CALL INPR(JPMAX,SP,TPS,DTP,D2P,JPMAX,SP,XS)
C
C
      WRITE(8)DTP,DTS,D2P,D2S
      IF(IBUG.EQ.1)WRITE(6,600)DTP,DTS
600   FORMAT(2(F10.5,5X))
100   CONTINUE
      REWIND 8
      RETURN
      END
      SUBROUTINE CALT(IWRITE,IDMP,IER)
C

C       THIS PROGRAM CALCULATES THE CAPTHETA TRANSFORMS
C       AT GIVEN GRID POINTS
C
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      CALL RDMID(NPTS,NL,NBLADE)
      CALL RGRID
C
C
C
C
      CALL INSEC(NPTS,NL,IER,IDMP)
      IF(IER.NE.0)GO TO 1000
C
C
10    CONTINUE
C
C
C
      CALL GETTHA(NPTS,NL,IDMP)
C
C
      CALL CALDR(IDMP)
C
C
      CALL CALDER(NBLADE,IDMP)
1000  CONTINUE
      RETURN
      END
      FUNCTION EXTR(TP,RP,R)
      REAL TP(50),RP(50)
      DELR=RP(2)-RP(1)
      DTDR=(TP(2)-TP(1))/DELR
      T= TP(1) + DTDR*(R-RP(1))
      EXTR=T
      RETURN
      END
      SUBROUTINE FILL(N,X,Y)
      REAL X(50),Y(50)
      COMMON /WQ/IWR
      COMMON/SCL/ SCALE
      READ(5,500)N,SCALE,ADJ
      IF(N.GT.0.AND.N.LE.50)GO TO 5
      WRITE(6,600)N
600   FORMAT(24H-POTENTIAL FILL ERROR N= ,I5)
5     CONTINUE
      IF(SCALE.LE.0.0)SCALE=1.0
      ADJ=ADJ/SCALE
      READ(5,501)(X(I),Y(I),I=1,N)
      WRITE(IWR,603)
      DO 10 I=1,N
      X(I)= X(I)/SCALE + ADJ
      Y(I)=Y(I)/SCALE
      WRITE(IWR,604)X(I),Y(I)
10    CONTINUE
500   FORMAT(I5,2F10.5)
501   FORMAT(2F10.5)
603   FORMAT(//'    SCALED RAW DATA INPUT   X,Y PAIRS')
604   FORMAT(10X,2F10.5)
      RETURN
      END
      SUBROUTINE GETTHA(NPTS,NL,IDMP)
C       THIS ROUTINE CALCULATES THE  THETA POSITIONS OF
C       SUCTION AND PRESSURE SURFACES AT THE INTERSECTION
C       POINTS
C
C
      COMMON /GRID/NX,NTH,NR,DQ,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      REAL MID
      COMMON /MIDSAV/MID(18,50,4)
      REAL ZI(50),RI(50),THP(50),THS(50),DT2SDS(50)
      REAL XIC(50),RIC(50),THPS(50),THSS(50),DTSDS(50),DT2PDS(50)
      REAL X(50),R(50),TP(50),TS(50),SP(50),SS(50),DTPDS(50)
      EQUIVALENCE (XIC,SS),(RIC,SP)
      NR1=NR-1
      IBUG=0
      IF(IDMP.LT.0  .OR. IDMP.EQ.3)IBUG=1
      DO 100 K=N1,N2
      READ(8) XIC,RIC
      READ(2)X,R
      DO 20 L=1,NL
      DO 10 J=1,NPTS
      ZI(J)=MID(L,J,1)
      RI(J)=MID(L,J,2)
      THP(J)=MID(L,J,3)
      THS(J)=MID(L,J,4)
10    CONTINUE
C
C
      XQ=XIC(L)
      CALL INTER(XQ,DX,D2X,NPTS,ZI,THP,T)
      THPS(L)=T
      CALL INTER(XQ,DX,D2X,NPTS,ZI,THS,T)
      THSS(L)=T
20    CONTINUE
      IF(IBUG.EQ.1)WRITE(6,610)K,NL,(THPS(L),THSS(L),L=1,NL)
610   FORMAT(2I5,2(F10.5,5X))
C
C
      DO 50 L=1,NR
      RQ = R(L)
C
      IF(RQ.GE.RIC(1))GO TO 25
C
C       USE LINEAR EXTRAPOLATION TO FILL OUT BLADE DATA
C
      WRITE(6,600)K,L
600   FORMAT(' EXTRAPOLATION NEEDED FOR J,L ',2I5)
      T=EXTR(THPS,RIC,RQ)
      TP(L)=T
      T=EXTR(THSS,RIC,RQ)
      TS(L)=T
      GO TO 50
C
C
25    CONTINUE
      CALL INTER(RQ,DX,D2X,NL,RIC,THPS,T)
      TP(L)=T
      CALL INTER(RQ,DX,D2X,NL,RIC,THSS,T)
      TS(L)=T
50    CONTINUE
C
C       CALCULATE ARC LENGTH ALONG CAPX=CONST. LINES
C
      CALL ARC(NR,X,R,TS,SS)
      CALL ARC(NR,X,R,TP,SP)
C
C       CALCULATE DTHETA/DS
C
      CALL INPR(NR,SS,TS,DTSDS,DT2SDS,NR,SS,THPS)
      CALL INPR(NR,SP,TP,DTPDS,DT2PDS,NR,SP,THPS)
C
C
      WRITE(4)TS,TP,DTPDS,DTSDS,DT2PDS,DT2SDS
C
      IF(IBUG.EQ.1)WRITE(6,620)TP,TS,DTPDS,DTSDS
620   FORMAT(4(F10.5,5X))
100   CONTINUE
C
C
      REWIND 4
      REWIND 2
      REWIND 8
      RETURN
      END
      SUBROUTINE INDP(X,DRDX,DR2,N,Z,R,RV)
      REAL R(50),Z(50)
      IF(X.GT.Z(1)) GO TO 10
      DRDX=0.
      DR2=0.
      RV=R(1)
      RETURN
10    CONTINUE
      IF(X.LT.R(N))GO TO 20
      DRDX=0.
      DR2=0.
      RV=R(N)
15    RETURN
20    CONTINUE
      DO 100 I=2,N
      I1=I-1
      IF(X.GT.Z(I))GO TO 100
      DRDX=(R(I)-R(I1))/(Z(I)-Z(I1))
      DR2=0.
      RV=R(I1)+DRDX*(X-Z(I1))
      GO TO 15
100   CONTINUE
      DRDX=0.
      DR2=0.
      RV=R(N)
      RETURN
      END
      SUBROUTINE INPR(N1,X,Y,D,D2,N2,Z,YI)
      REAL X(50),Y(50),D(50),D2(50),Z(50),YI(50)
      REAL S(50),EM(50),WORK(2,50)
      NA=N1
      NB=N2
      CALL SPLINT(NA,NB,X,Y,Z,YI,D,D2,S,EM,WORK)
C CALCULATE THE SECOND DERIVATIVES WITH SECOND ORDER ACCURACY
      NA1=NA-1
      DO 10 I=2,NA1
      DIP1=X(I+1)-X(I)
      DIM1=X(I)-X(I-1)
      D2(I)=D(I)*(DIP1-DIM1)/DIP1/DIM1
      PRE=D(I+1)*DIM1/DIP1-D(I-1)*DIP1/DIM1
      D2(I)=D2(I)+PRE/(DIP1+DIM1)
10    CONTINUE
      D2(1)=0.0
      D2(NA)=0.0
      RETURN
      END
      SUBROUTINE INSEC(NPTS,NL,IER,IDMP)
      COMMON /MIDSAV/ MID(18,50,4)
      REAL MID
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON /GRID/NX,NTH,NP,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      REAL ZI(50),RI(50),TI(50)
      REAL TQ(50)
      REAL X(50),R(50),XIC(50),RIC(50)
      TOL=0.0001
      NRSAV=NP
      IBUG=0
      IF(IDMP.LT.0  .OR. IDMP.EQ.2)IBUG=1
      IER=0
      NR=NRSAV+2
      NR1=NR-1
      DO 100 J=N1,N2
      JP= J-N1+1
      READ(2)X,R,TQ,TQ,TQ,TQ
C
C              EXTEND X,R TO FORCE INTER SECTIONS
C
      R(NR)=1.3
      X(NR)=X(NRSAV)
      IR=NR-1
4     IR1=IR-1
      X(IR)=X(IR1)
      R(IR)=R(IR1)
      IR=IR-1
      IF(IR.GE.2)GO TO 4
      R(1)=0.
      X(1)=X(2)
C
C
      IT=0
      DO 95 KP=1,NL
      DO  3  I=1,NPTS
      ZI(I)=MID(KP,I,1)
      RI(I)=MID(KP,I,2)
3     CONTINUE
      ICOUNT=0
C
C       TEST FOR NO POSSIBLE INTERSECTIONS
C
      R1=RI(1)
      IF(R1.GT.R(1).AND. R1.LT.R(NR))GO TO 5
      WRITE(6,600)J,KP
600   FORMAT('  NO POSSIBLE INTERSECTION FOR J,KP ',2I5)
      IER=1
      GO TO 200
C
C
5     CONTINUE
C
C       GUESS X INTER. POINT
C
      XQ=X(1)
C
C       WHAT R POSITION ON BLADE LINE DOES THIS CORRESPOND TO
C
10    CALL INTER(XQ,DX,D2X,NPTS,ZI,RI,RQ)
C
C       WHAT X ON GRID INPUT LINE => RQ
C
      CALL INTER(RQ,DX,D2X,NR,R,X,XNEW)
C
C              FORCE AT LEAST 1 ITERATION
C
      IF(ICOUNT.EQ.0)XQ=-4000.
      IF(ABS(XNEW-XQ).LE.TOL)GO TO 90
C
C       NO MATCH YET
C
      ICOUNT=ICOUNT+1
      IF(ICOUNT.LE.10)GO TO 20
      WRITE(6,601)J,KP
601   FORMAT('  NO MATCH FOUND FOR J,KP ',2I5)
      IER=1
      GO TO 200
20    CONTINUE
C
C              GUESS NEW X
      XQ=XNEW
      GO TO 10
C
C       INTERSECTION FOUND
C
90    CONTINUE
      IT=IT+1
      XIC(IT)=XQ
      RIC(IT)=RQ
95    CONTINUE
      WRITE(8)XIC,RIC
      IF(IBUG.EQ.1)WRITE(6,605)J,(XIC(I),RIC(I),I=1,IT)
605   FORMAT(I5,/,2(F10.5,5X))
100   CONTINUE
200   CONTINUE
      REWIND 2
      REWIND 8
      RETURN
      END
      SUBROUTINE INTER(X,DYDX,D2YDX,NPT,XH,RH,Y)
      REAL XH(50),RH(50),EM(50),S(50),WORK(2,50),TP(2)
      REAL D(2),D2(2),R(2)
      DIMENSION XE(2),YE(2)
C.... NEW VERSION OF INTER, WILL EXTRAPOLATE
C
      IF(X.LE.XH(NPT))GO TO 10
C.... X OUT OF RANGE TO RIGHT.
      NP1 = NPT-1
      XE(1) = XH(NP1)
      XE(2) = XH(NPT)
      YE(1) = RH(NP1)
      YE(2) = RH(NPT)
      GO TO 20
10    IF(X.GE.XH(1)) GO TO 30
C.... X OUT OF RANGE TO LEFT.
      XE(1) = XH(2)
      XE(2) = XH(1)
      YE(1) = RH(2)
      YE(2) = RH(1)
20    Y = EXTR(YE,XE,X)
      DYDX=0.
      D2YDX=0.
      RETURN
C
30    CONTINUE
      NP=NPT
      NI=1
      TP(1)=X
      CALL SPLINT(NP,NI,XH,RH,TP,R,D,D2,S,EM,WORK)
      Y=R(1)
      DYDX=D(1)
      D2YDX=D2(1)
      RETURN
      END
      SUBROUTINE RDMID(NPTS,NL,NBLADE)
C
C       THIS SUBROUTINE READS THE BLADE GEOMETRY DATA
C
      REAL MID
      COMMON /MIDSAV/MID(18,50,4)
      COMMON/FILEWR/IW
      REAL TITLE(90)
      REAL XX(50),Y(50),Z(50)
      PI=3.1416
      READ(5,510)(TITLE(I),I=1,80)
      READ(5,520)NPTS,NL,SCALE,IM,NBLADE,ADJ
      BIGTHA=2.*PI/NBLADE
      IF(SCALE.LE.0.0) SCALE = 1.0
C
C       READ Z,R,THETA OF MIDCHORD LINE AND TANG. THICKNESS
C
      DO 150 K=1,4
      DO 150 I=1,NL
      READ(5,627)(MID(I,J,K),J=1,NPTS)
627   FORMAT(8F10.3)
150   CONTINUE
      IF(IM.EQ.1) GO TO 200
C
C       CALCULATE THETA PRESSURE SUCTION SURFACES
C
      DO 170 L=1,NL
      DO 165 J=1,NPTS
      RTHMID=MID(L,J,3)*MID(L,J,2)
C
C       CORRECT THETA DIRECTION DIFFERENCE BETWEEN NASA AND MIT
C
      TH=-MID(L,J,3)
      TT=MID(L,J,4)/2.
      MID(L,J,3)=TH+TT
      MID(L,J,4)=TH-TT
165   CONTINUE
170   CONTINUE
C
C
200   CONTINUE
C
C       OPTION 2 OF INPUT SHOULD BE TAKEN CARE OF ALREADY I.E.
C       PREVIOUS CODE FOR OPTION 1 SHOULD CREATE OPTION 2 INPUT
C
C       SCALE R AND Z
C       ADD 2*PI/NBLADE TO SUCTION SURFACE THETA
C
      TZ=ADJ/SCALE
      DO 240 L=1,NL
      DO 230 J=1,NPTS
      DO 220 I=1,2
      MID(L,J,I)=MID(L,J,I)/SCALE
220   CONTINUE
      MID(L,J,4)=MID(L,J,4)+BIGTHA
      MID(L,J,1)=MID(L,J,1)+TZ
230   CONTINUE
240   CONTINUE
C
C       DISPLAY BLADE DATA

C
      IF(IM.EQ.0)WRITE(IW,515)(TITLE(I),I=1,80)
      DO 330 L=1,NL
      WRITE(IW,601)L
      DO 320 J=1,NPTS
      RAD= MID(L,J,2)
      RTHPS=RAD*MID(L,J,3)
      RTHSS=RAD*MID(L,J,4)
      RTHSS=RTHSS-BIGTHA*RAD
      Y(J)=RTHPS
      Z(J)=RTHSS
      X=MID(L,J,1)
      XX(J)=X
      WRITE(IW,600)L,J,X,RAD,RTHPS,RTHSS
320   CONTINUE
      WRITE(IW,602)
      CALL SMOOTH(XX,Y,NPTS,IW)
      WRITE(IW,603)
      CALL SMOOTH(XX,Z,NPTS,IW)
330   CONTINUE
510   FORMAT(80A1)
515   FORMAT(///,1X,80A1)
520   FORMAT(2I5,F10.5,2I5,F10.5)
600   FORMAT(2I5,4F10.5)
601   FORMAT(/,' INPUT DATA FOR SECTION ',I5,/,
     1T5,'L',T10,'J',T20,'X',T30,'R',T36,'RTHPS',T46,'RTHSS')
602   FORMAT(/' SECOND DERIVATIVE CHECK FOR PRESSURE SURFACE')
603   FORMAT(/' SECOND DERIVATIVE CHECK FOR SUCTION SURFACE')
      RETURN
      END
      SUBROUTINE RECLUS(KX,NR,X,R)
      REAL KX,X(50),R(50),NEWR(50),NEWX(50)
      R1=R(1)
      DELR=2.0/(NR-1)
      D=R(NR)-R(1)
      DO 100 L=1,NR
      T=-1.0 + DELR*(L-1)
      CALL STRTCH(T,T1,KX,AX,BX)
      T2=AX*SINH(KX*T)+BX*ATAN(SINH(KX*T))
      T3=D*(T2+0.5)+R1
      NEWR(L)=T3
100   CONTINUE
      DO 200 L=1,NR
      R1=NEWR(L)
      CALL INTER(R1,T1,T2,NR,R,X,X1)
      NEWX(L)=X1
200   CONTINUE
      DO 300 L=1,NR
      R(L)=NEWR(L)
      X(L)=NEWX(L)
300   CONTINUE
      RETURN
      END
      SUBROUTINE RGRID
C
C       THIS ROUTINE READS AND STORES THE POSITION OF CONSTANT
C       CAP-X LINES IN PHYSICAL SPACE
C
      REAL X(50),R(50),XX(50),XR(50),RR(50),RX(50)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      DO 100 J=1,NX
      DO 10 L=1,NR
      READ(3)X(L),R(L),RR(L),RX(L),XX(L),XR(L)
10    CONTINUE
C
      READ(3)T1,T2
      READ(3)T1,T2
      IF(J.LT.N1)GO TO 100
      IF(J.GT.N2)GO TO 200
      WRITE(2)X,R,XX,XR,RR,RX
100   CONTINUE
200   CONTINUE
      REWIND 2
      REWIND 3
      RETURN
      END
      SUBROUTINE SECDER(X1,X2,X3,D1,D2,D3,DSEC)
C
C       THIS SUBROUTINE CALCULATES THE SECOND DERIVATIVES
C       FOR THE HUB AND TIP SLOPES WITH SECOND ORDER ACCU-
C       RACY FROM THE DERIVATIVES FROM SPLINT
C
      DEL1=X2-X1
      DEL3=X3-X2
      DSEC=D2*(DEL3-DEL1)/DEL1/DEL3-DSEC
      RETURN
      END
      SUBROUTINE SETR(DRHDX,DRTDX,RV,TH,TT,D2H,D2T)
      REAL RV(4)
      PHI=ATAN(DRHDX)
      RV(1)=COS(PHI)
      RV(2)=-SIN(PHI)
      PHI=ATAN(DRTDX)
      RV(3)=-COS(PHI)
      RV(4)=SIN(PHI)
      D2H=-D2H
      TH=(1.+DRHDX**2)**1.5
      TH=D2H/TH
      TT=(1.+DRTDX**2)**1.5
      TT=D2T/TT
      RETURN
      END
      SUBROUTINE SMOOTH(X,Y,M,IWR)
      REAL X(50),Y(50),YI(50),D(50),D2(50),S(50),EM(50),WORK(2,50)
      N=M
C       THIS ROUTINE COMPUTES THE SECOND DERIVATIVES OF THE INPUT
C       DATA AS AN AID IN DETERMINING INPUT ERRORS
      CALL SPLINT(N,N,X,Y,X,YI,D,D2,S,EM,WORK)
C
C
      WRITE(IWR,500)
500   FORMAT(/,1X,T5,'J',T15,'X',T25,'Y',T32,'DYDX     D2YDX')
      DO 50 J=1,M
      WRITE(IWR,501)J,X(J),Y(J),D(J),D2(J)
50    CONTINUE
501   FORMAT(I5,5F10.4)
      RETURN
      END
      SUBROUTINE SPLINT(N, MAX, X, Y, Z, YINT, DYDX, D2YDX, S, EM, WORK)
C
C        CALLING SEQUENCE
C              N - THE NUMBER OF SPLINE POINTS (MINIMUM 2)
C              MAX - THE NUMBER OF INTERPOLATION POINTS
C              X - VECTOR OF X COORDINATES OF N SPLINE POINTS
C              Y - VECTOR OF Y COORDINATES OF N SPLINE POINTS
C              Z - VECTOR OF X-COORDINATES OF INTERPOLATION POINTS
C              YINT - VECTOR OF Y COORDINATES CORRESPONDING TO POINTS Z
C              DYDX - VECTOR OF DY/DX         CORRESPONDING TO POINTS Z
C              D2YDX- VECTOR OF D2Y/DX2       CORRESPONDING TO POINTS Z
C              S,EM - ARRAYS OF DIMENSION N CALCULATES BY SPLINT
C              WORK - ARRAY USED BY SPLINT AS WORKING SPACE
C
      DIMENSION X(N), Y(N)
      DIMENSION Z(MAX), YINT(MAX), DYDX(MAX), D2YDX(MAX)
      DIMENSION S(N), EM(N)
      DIMENSION WORK(2,N)
1000  FORMAT(17H0OUT OF RANGE X =, G14.6)
1001  FORMAT('0   ERROR: SPLINT ENTERED WITH N LESS THAN TWO')
C***********************************************************************
      IF(N.LT.2) GO TO 900
      DO 10 I=2,N
10    S(I)=X(I)-X(I-1)
      NO=N-1
      WORK(1,1) =-.5
      WORK(2,1)=0.
      IF(2.GT.NO) GO TO 35
      DO 30 I=2,NO
      SI=S(I)
      SI1=S(I+1)
      F=(Y(I+1)-Y(I))/SI1 - (Y(I)-Y(I-1))/SI
      SID6=0.16666667*SI
      W=(SI+SI1)*.3333333-SID6*WORK(1,I-1)
      WORK(1,I)=(SI1*.1666667)/W
      WORK(2,I)=(F-SID6*WORK(2,I-1))/W
30    CONTINUE
35    EM(N)=(.5*WORK(2,NO))/(1.+.5*WORK(1,NO))
      DO 40 I=2,N
      K=N+1-I
40    EM(K)=WORK(2,K)-WORK(1,K)*EM(K+1)
      IF(MAX.LT.1) RETURN
      DO 150 I=1,MAX
      ZI=Z(I)
      K=2
      IF(ZI-X(1)) 50, 100, 55
50    IF(ZI.LT.(1.1*X(1)-.1*X(2))) WRITE(6,1000) ZI
      GO TO 100
55    K=N
      IF(ZI-X(N)) 65, 100, 60
60    IF(ZI.GT.(1.1*X(N)-.1*X(N-1))) WRITE(6,1000) ZI
      GO TO 100
65    MN1=2
      MX=N
70    K=(MX+MN1)/2
      IF(MX.EQ.MN1) GO TO 100
      IF(ZI-X(K)) 75, 100, 80
75    MX=K
      GO TO 70
80    MN1=K+1
      GO TO 70
100   XZ=X(K)-ZI
      XZ2=XZ*XZ
      ZX=ZI-X(K-1)
      ZX2=ZX*ZX
      SK=S(K)
      SK2=SK*2.
      SKD6=SK*.16666667
      EMK=EM(K)
      EMK1=EM(K-1)
      YKDS=Y(K)/SK
      YK1DS=Y(K-1)/SK
      YINT(I)=(EMK1*XZ2*XZ+EMK*ZX2*ZX)*.1666667/SK + (YKDS-EMK*SKD6)*ZX
     *   + (YK1DS-EMK1*SKD6)*XZ
      DYDX(I)= (EMK*ZX2-EMK1*XZ2)/SK2 + YKDS - YK1DS - (EMK-EMK1)*SKD6
      D2YDX(I)=(EMK1*XZ+EMK*ZX)/SK
150   CONTINUE
      RETURN
900   WRITE(6,1001)
      RETURN
      END
      SUBROUTINE STRTCH(X,DXDX,KX,AX,BX)
      REAL KX
      SH=SINH(KX)
      CH=COSH(KX)
      ASTAR=1./CH/2./KX
      BSTAR=SH*TANH(KX)/2./KX
      AA=ASTAR*SH+BSTAR*ATAN(SH)
      AA=2.*AA
      AX=ASTAR/AA
      BX=BSTAR/AA
      SH=SINH(KX*X)
      CH=COSH(KX*X)
      T=1.+SH**2
      T=AX+BX/T
      DXDX=1./KX/CH/T
      RETURN
      END
      SUBROUTINE SVGRID
C.... SVGRID SAVES ALL THE INFORMATION NEEDED FOR THE GRID AND SOLUTION
C     PLOT PACKAGES IN A FORMATTED FILE ON UNIT 12.
      COMMON/GRID/  NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/ N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/HUBTIP/XH(50),RH(50),XT(50),RT(50),NRH,NRT
      COMMON/LEAD/  ZLE(30),RLE(30),ZTE(30),RTE(30),NZLE,NZTE,A1
      COMMON/DAMP/  IDAM,ZL(20),RL(20),ZU(20),RU(20),NZL,NZU,C1,C2,C3
      COMMON/ROTOR/ TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON/SCL/   SCALE
      DIMENSION  X(100,18),  R(100,18)
      DIMENSION TS(100,18), TP(100,18)
      DIMENSION TSR(50),TPR(50)
      DIMENSION DTSDZR(50),DTPDZR(50)
      DIMENSION DTDZ(20,2)
      DIMENSION DU2(50)
      PITCH = 6.283185/FLOAT(NBLADE)
C.... WRITE GENERAL INFORMATION.
      IF(IEXPLE.GT.0) GO TO 10
      NZL = 0
      NZU = 0
10    CONTINUE
      WRITE(12)NBLADE,NX,NTH,NR,N1,N2,KUTTA,IEXPLE,IEXPTE,
     1NRH,NRT,NZLE,NZTE,NZL,NZU
      WRITE(12)SCALE
C.... WRITE HUB, TIP, LE, TE, AND DAMPER COORDINATES.
      WRITE(12)( XH(I),I=1,NRH )
      WRITE(12)( RH(I),I=1,NRH )
      WRITE(12)( XT(I),I=1,NRT )
      WRITE(12)( RT(I),I=1,NRT )
      WRITE(12)(ZLE(I),I=1,NZLE)
      WRITE(12)(RLE(I),I=1,NZLE)
      WRITE(12)(ZTE(I),I=1,NZTE)
      WRITE(12)(RTE(I),I=1,NZTE)
      IF(IEXPLE.LE.0) GO TO 100
      WRITE(12)( ZL(I),I=1,NZL )
      WRITE(12)( RL(I),I=1,NZL )
      WRITE(12)( ZU(I),I=1,NZU )
      WRITE(12)( RU(I),I=1,NZU )
100   CONTINUE
C.... READ THE X-R COORDINATES.
      DO 410 J=1,NX
      DO 400 L=1,NR
400   READ(3) X(J,L),R(J,L),DUM,DUM,DUM,DUM
C.... SKIP 2 RECORDS, OR 3 IF A DAMPER IS PRESENT.
      READ(3) DUM
      READ(3) DUM
      IF(IEXPLE.LE.0) GO TO 410
      READ(3) DUM
410   CONTINUE
      REWIND 3
C.... READ THE THETA COORDINATES.
      NR1 = NR-1
      DO 420 J=N1,N2
      READ(9) DU2,DU2,TSR,TPR,DTPDZR,DTSDZR,DU2,DU2
      READ(9) DUM,DUM,DUM
      READ(9) DUM,DUM
      DO 420 L=1,NR
      TS(J,L) = TSR(L)
      TP(J,L) = TPR(L)
C.... CALCULATE THE MEAN CAMBER SLOPES.
C.... DTDZ(L,1)=>L.E., DTDZ(L,2)=>T.E.
      IF(J.NE.N1.AND.J.NE.N2) GO TO 420
      DTHDZ = 0.5*(DTSDZR(L)+DTPDZR(L))
      IF(J.EQ.N1) DTDZ(L,1) = DTHDZ
      IF(J.EQ.N2) DTDZ(L,2) = DTHDZ
420   CONTINUE
      REWIND 9
C.... CALCULATE UP AND DOWNSTREAM THETAS BY LINEAR EXTRAPOLATION.
      N1M1 = N1-1
      N2P1 = N2+1
      DO 450 L=1,NR
      XLE  =  X(N1,L)
      XTE  =  X(N2,L)
      TSLE = 0.5*(TS(N1,L)+TP(N1,L)+PITCH)
      TSTE = 0.5*(TS(N2,L)+TP(N2,L)+PITCH)
      DLEDZ= DTDZ(L,1)
      DTEDZ= DTDZ(L,2)
      DO 430 J=1,N1M1
      TS(J,L) = DLEDZ*(X(J,L)-XLE)+TSLE
430   TP(J,L) = TS(J,L)-PITCH
      DO 440 J=N2P1,NX
      TS(J,L) = DTEDZ*(X(J,L)-XTE)+TSTE
440   TP(J,L) = TS(J,L)-PITCH
450   CONTINUE
C.... WRITE X, R, TS, AND TP.
      WRITE(12)(( X(J,L),J=1,NX),L=1,NR)
      WRITE(12)(( R(J,L),J=1,NX),L=1,NR)
      WRITE(12)((TS(J,L),J=1,NX),L=1,NR)
      WRITE(12)((TP(J,L),J=1,NX),L=1,NR)
      RETURN
      END
      SUBROUTINE TRNR(X,R,RR,RZ,CAPR,RHUB)
      COMMON/HUBTIP/XH(50),RH(50),XT(50),RT(50),NRH,NRT
C
      CALL INTER(X,DRHDX,D2RHDX,NRH,XH,RH,RHUB)
      CALL INTER(X,DRTDX,D2RTDX,NRT,XT,RT,RTIP)
C
      DELR=RTIP-RHUB
      CAPR=(R-RHUB)/DELR
      RR=1./DELR
      T1=DELR*DRHDX
      T2=(R-RHUB)*(DRTDX-DRHDX)
      T=DELR**2
      RZ=-1.*(T1+T2)/T
      RETURN
      END
      SUBROUTINE TRNX(X,R,ZZ,ZR,CAPX,ZL)
      COMMON/LEAD/ZLE(30),RLE(30),ZTE(30),RTE(30),NZLE,NZTE,A1
      COMMON/PLACE/JX
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
C
      CALL INTER(R,DZLE,T,NZLE,RLE,ZLE,ZL)
      CALL INTER(R,DZTE,T,NZTE,RTE,ZTE,ZT)
C
      DELZ=ZT-ZL
      ZZ=1./DELZ
      T1=DELZ*DZLE
      T2=(X-ZL)*(DZTE-DZLE)
      T3=DELZ**2
      ZR=-1.*(T1+T2)/T3
      CAPX=(X-ZL)/DELZ
      CAPX=CAPX-0.5
      IF(JX.LT.N1)GO TO 10
      IF(JX.GT.N2)GO TO 20
      RETURN
10    CONTINUE
      ZSTAR=X-ZL
      ZSTARP=+DZLE
      SIXP=1.
      GO TO 30
20    CONTINUE
      ZSTAR=ZT-X
      ZSTARP=-DZTE
      SIXP=-1.
30    CONTINUE
      ZL0=ZLE(1)
      ZT0=ZTE(1)
      ZLP=ZL0-ZL
      ZTP=ZT-ZT0
      EX=EXP(A1*ZSTAR)
      DEX=A1*ZSTARP*EX
      TOP=(X-ZL0)+ZLP*EX
      AA1=ZL0-ZL
      AA2=ZT0 - ZL0
      AA3=ZT - ZT0 + ZL0 - ZL
      BOT=AA2+AA3*EX
      B1=X-ZL0
      TOP=B1+AA1*EX
      CAPX=TOP/BOT - 0.5
      T1=(AA3+SIXP*AA1*AA2*A1-SIXP*
     #     AA3*A1*(X-ZL0))*EX
      ZZ=(AA2+T1)/(BOT*BOT)
      T1=DZLE*(TOP-BOT)-DZTE*TOP
      T2=ZSTARP*A1*(AA3*B1-AA1*AA2)
      ZR=EX*(T1+T2)/(BOT*BOT)
      RETURN
      END
      SUBROUTINE XRGRID(IWRITE,IDMP,IER,IAN)
C
C
C       THIS PROGRAM GENERATES THE X-R GRID.
C       HAS BEEN MODIFIED IN ORDER TO COMPUTE DERIVATIVES AND
C       ALSO SOME OTHER MINOR MODIFICATIONS
C
C              CHECKED 10-MAY-79   **GHH**
C
C
      REAL KX
      REAL ZZ(100),Z(100)
      REAL DXDZZ(100)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/DAMP/IDAM,ZLL(20),RL(20),ZU(20),RU(20),NZL,NZU,C1,C2,C3
      COMMON /WQ/IWR
      COMMON /HUBTIP/XH(50),RH(50),XT(50),RT(50),NRH,NRT
      COMMON /LEAD/ZLE(30),RLE(30),ZTE(30),RTE(30),NZLE,NZTE,A1
      COMMON/PLACE/JX
      REAL PH(50),PR(50),RRS(50),RXS(50),ZXS(50),ZRS(50)
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON /XRTSAV/XSAVE(50,100,3),RSAVE(50,100,3)
      COMMON/TTSAVE/TTSAV(50,50,3),TRSAV(50,50),TXSAV(50,50),
     *   TPSAV(50,100),TSSAV(50,100)
      COMMON /MESHSP/SPHALF,SKR,SKT
      COMMON/DAMPPS/IDLOW,IDHIGH
C
C
      NXP=NX
      NR1=NR-1
      NR2=NR-2
      IWR=IWRITE
      NO=(N1+N2)/2
C
C       GET INPUT DATA
C
      READ(5,503)KX,A1,TOL,IAB
503   FORMAT(3F10.5,I5)
      WRITE(IWR,605)
605   FORMAT(/,' HUB, TIP, L.E., AND T.E. COORDINATES.')
      CALL FILL(NA,XH,RH)
      NPP=NA
      NRH=NA
      CALL FILL(NA,XT,RT)
      NRT=NA
      CALL FILL(NZLE,ZLE,RLE)
      CALL FILL(NZTE,ZTE,RTE)
C
C
      IDAM=0
      IF(IEXPLE.NE.0)IDAM=1
      IF(IDAM.EQ.0)GO TO 15
      WRITE(IWR,606)
606   FORMAT(/,' LOWER AND UPPER DAMPER INPUT COORDINATES.')
      CALL FILL(NZL,ZLL,RL)
      CALL FILL(NZU,ZU,RU)
      READ(5,503)C1
      C2 = C1*(NR-2)
      IDLOW = INT(C2)+1
      IDHIGH = IDLOW+1
      WRITE(IWR,610) IDLOW,IDHIGH
610   FORMAT(//'DAMPER IS BETWEEN RADIAL GRID PLANES NO.',I5,'AND',I5)
C
C
C
C
C
15    CONTINUE
C
      WRITE(IWR,607)
607   FORMAT(/,' HUB, TIP, L.E., AND T.E. DERIVITIVES.')
C
      CALL SMOOTH(XT,RT,NRT,IWR)
      CALL SMOOTH(XH,RH,NRH,IWR)
      CALL SMOOTH(RLE,ZLE,NZLE,IWR)
      CALL SMOOTH(RTE,ZTE,NZTE,IWR)
      IF(IDAM.EQ.0)GO TO 5
      WRITE(IWR,608)
608   FORMAT(/,' LOWER AND UPPER DAMPER DERIVITIVES.')
      CALL SMOOTH(ZLL,RL,NZL,IWR)
      CALL SMOOTH(ZU,RU,NZU,IWR)
C
C
5     CONTINUE
C       CALCULATE   X TO ZZ TRANSFORM
C
      A=0.
      B=0.
      DO 10 I=1,NX
      X= DX*(I-NO)
      CALL STRTCH(X,DXDY,KX,AX,BX)
      Z1=SINH(KX*X)
      ZZ(I)=Z1*AX +BX*ATAN(Z1)
      DXDZZ(I)=DXDY
10    CONTINUE
C
C       COMPUTE XLAST AND RHUB LAST FOR FIRST GRID POINT
C
      RHUBL=RH(1)
      CALL INTER(RHUBL,T1,T2,NZLE,RLE,ZLE,ZL)
      CALL INTER(RHUBL,T1,T2,NZTE,RTE,ZTE,ZT)
      CAPX=ZZ(1)
      XLAST=(CAPX + 0.5)*(ZT-ZL) + ZL
      WRITE(IWR,609)
609   FORMAT(/,' X-R GRID COORDINATES AND METRICS.')
C
      WRITE(IWR,600)KX
600   FORMAT(' THE AXIAL PACKING FACTOR IS ',F6.3)
C
C       ITERATION ALONG X STARTS HERE
C
      DO 55 J=1,NX

601   FORMAT(//,T5,'AXIAL',T15,'RADIAL',T25,'AXIAL',T35,
     * 'RADIAL',/,T5,'STATION',T15,'STATION',T25,'POSITION',
     * T35,'POSITION',T48,'XX',T58,'XR',T68,'RR',T78,'RX')
602   FORMAT(T7,I5,T17,I5,T25,F8.4,T35,F8.4,T45,F7.3,T55,F7.3,
     * T65,F7.3,T75,F7.3)
      JX=J
      CAPX=ZZ(J)
C
C       GUESS X AND R
C
      XQ=XLAST
      RQ=RHUBL
C
C
      DO 50 L=1,NR
      IT=0
      CAPR=DR*(L-1)
      CAPR=CAPR-DR*SPHALF
C
C       COMPUTE NEW CAPR AND CAPX
C
20    CALL TRNR(XQ,RQ,RR,RZ,CAPRQ,RHUB)
      CALL TRNX(XQ,RQ,ZX,ZR,CAPXQ,ZL)
C
C       CONVERGED _
C
      IT=IT+1
      IF(IT.LT.25)GO TO 35
      WRITE(6,650)JX
650   FORMAT(//'  ITERATION NOT CONVEREGED AT ',I5)
      IF(IAB.EQ.1)GO TO 45
      IF(IAB.EQ.0)GO TO 500
C
C
      IT=0
35    CONTINUE
      TT1=CAPX-CAPXQ
      DCAP=CAPX
      IF(CAPX.LT.1E-10)DCAP=1.0
      TT1=ABS(TT1/DCAP)
      TT2=CAPR-CAPRQ
      DCAP=CAPR
      IF(CAPR.LT.1E-10)DCAP=1.0
      TT2=ABS(TT2/DCAP)
      IF(IT.EQ.1)GO TO 40
      IF(TT1.GE.TOL)GO TO 40
      IF(TT2.LE.TOL)GO TO 45
C
C       NOT CONVERGED; CALCULATE NEW XQ,RQ
C
40    CONTINUE
      XQS=XQ
      RQS=RQ
C       XQ =(CAPX+0.5)/ZX +ZL
      XQ=(CAPX-CAPXQ)/ZX + XQ
      RQ=(CAPR-CAPRQ)-RZ*(XQ-XQS)
      RQ=RQ/RR
      RQ=RQ+RQS
      GO TO 20
C
C       CONVERGED
C
45    CONTINUE
      IF(L.NE.2) GO TO 46
      XLAST= XQ
      RHUBL= RQ
46    PH(L)= XQ
      PR(L)= RQ
      RQ=RQ+DR/RR
50    CONTINUE
      IF(SKR.NE.0)CALL RECLUS(SKR,NR,PH,PR)
      DO 55 L=1,NR
      XSAVE(L,J,1)=PH(L)
      RSAVE(L,J,1)=PR(L)
55    CONTINUE
      CALL XXM(IAN)
      DO 100 J=1,NX
      WRITE(IWR,601)
C
C       CALCULATE INTERSECTIONS, ANGLES AND CURVATURES
C
      DO 60 L=1,NR
C
      XQ=XSAVE(L,J,1)
      RQ=RSAVE(L,J,1)
      PH(L)=XQ
      PR(L)=RQ
      RR=RSAVE(L,J,2)
      RX=RSAVE(L,J,3)
      ZR=XSAVE(L,J,2)
      ZX=XSAVE(L,J,3)
      IF(IAN.EQ.1)ZX=ZX*DXDZZ(J)
      IF(IAN.EQ.1)ZR=ZR*DXDZZ(J)
      RRS(L)=RR
C
C       SAVE AXIAL AND RADIAL POSITIONS
C
      RXS(L)=RX
      ZXS(L)=ZX
      ZRS(L)=ZR
      WRITE(IWR,602)J,L,XQ,RQ,ZX,ZR,RR,RX
60    CONTINUE
C
      XHP=PH(1)+(PH(2)-PH(1))*SPHALF
      XTP=PH(NR)-(PH(NR)-PH(NR-1))*SPHALF
C
C       CALCULATE FIRST DERIVATIVE OF HUB AND TIP BY SPLINT
C
      CALL INTER(XTP,DRTDX,D2RTDX,NRT,XT,RT,RTIP)
      CALL INTER(XHP,DRHDX,D2RHDX,NRH,XH,RH,RHUB)
C
C       CALCULATE SECOND DERS. WITH SECOND ORDER ACCY.
C
      IF(J.EQ.1)GO TO 70
C
      DRH1=DRH2
      DRT1=DRT2
      XH1=XH2
      XT1=XT2
C
      DRH2=DRH3
      DRT2=DRT3
      XH2=XH3
      XT2=XT3
C
      DRH3=DRHDX
      DRT3=DRTDX
      XH3=XHP
      XT3=XTP
C
      IF(J.EQ.2)GO TO 70
      CALL SECDER(XH1,XH2,XH3,DRH1,DRH2,DRH3,D2RHDX)
      CALL SECDER(XT1,XT2,XT3,DRT1,DRT2,DRT3,D2RTDX)
      WRITE(3)DRH2,D2RHDX
      WRITE(3)DRT2,D2RTDX
C
C       FILL TRN.DAT
C
70    CONTINUE
      DO 75 L=1,NR
      XQ=PH(L)
      RQ=PR(L)
      RR=RRS(L)
      RX=RXS(L)
      ZX=ZXS(L)
      ZR=ZRS(L)
      WRITE(3)XQ,RQ,RR,RX,ZX,ZR
75    CONTINUE
C
      IF(J.NE.1.AND.J.NE.NX)GO TO 80
      WRITE(3)DRHDX,D2RHDX
      WRITE(3)DRTDX,D2RTDX
      DRH3=DRHDX
      DRT3=DRTDX
      XH3=XHP
      XT3=XTP
80    CONTINUE
C
C
C              COMPUTE DAMPER SLOPES
C
      IF(IDAM.EQ.0)GO TO 100
      IF(J.LT.IEXPLE) GO TO 100
      IF(J.GT.IEXPTE)GO TO 100
      XD=(PH(11)+PH(12))/2.
      CALL INDP(XD,DRDAML,D2DAML,NZL,ZLL,RL,R3)
      CALL INDP(XD,DRDAMU,D2DAMU,NZU,ZU,RU,R3)
      WRITE(3)DRDAML,D2DAML,DRDAMU,D2DAMU

100   CONTINUE
      REWIND 3
      RETURN
C
C
C
500   CONTINUE
      RETURN
      END
      SUBROUTINE XXM(IAN)
      COMMON /XRTSAV/XSAVE(50,100,3),RSAVE(50,100,3)
      COMMON/TTSAVE/TTSAV(50,50,3),TRSAV(50,50),TXSAV(50,50),
     *   TPSAV(50,100),TSSAV(50,100)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/PLACE/JX
      DX1=0.5/DX
      DR1=0.5/DR
      DO 100 J=1,NX
      JP=J+1
      JM=J-1
      JX=J
      DO 50 L=1,NR
      LP=L+1
      LM=L-1
      IF(IAN.EQ.1)GO TO 45
      IF(J.EQ.1)GO TO 10
      IF(J.EQ.NX)GO TO 15
      XJ=(XSAVE(L,JP,1)-XSAVE(L,JM,1))*DX1
      RJ=(RSAVE(L,JP,1)-RSAVE(L,JM,1))*DX1
      GO TO 20
10    XJ=(-3.*XSAVE(L,J,1)+4.*XSAVE(L,J+1,1)-XSAVE(L,J+2,1))*DX1
      RJ=(-3.*RSAVE(L,J,1)+4.*RSAVE(L,J+1,1)-RSAVE(L,J+2,1))*DX1
      GO TO 20
15    XJ=( 3.*XSAVE(L,J,1)-4.*XSAVE(L,J-1,1)+XSAVE(L,J-2,1))*DX1
      RJ=( 3.*RSAVE(L,J,1)-4.*RSAVE(L,J-1,1)+RSAVE(L,J-2,1))*DX1
20    IF(L.EQ.1) GO TO 30
      IF(L.EQ.NR)GO TO 35
      XL=(XSAVE(LP,J,1)-XSAVE(LM,J,1))*DR1
      RL=(RSAVE(LP,J,1)-RSAVE(LM,J,1))*DR1
      GO TO 40
30    XL=(-3.*XSAVE(L,J,1)+4.*XSAVE(L+1,J,1)-XSAVE(L+2,J,1))*DR1
      RL=(-3.*RSAVE(L,J,1)+4.*RSAVE(L+1,J,1)-RSAVE(L+2,J,1))*DR1
      GO TO 40
35    XL=(3.*XSAVE(L,J,1)-4.*XSAVE(L-1,J,1)+XSAVE(L-2,J,1))*DR1
      RL=(3.*RSAVE(L,J,1)-4.*RSAVE(L-1,J,1)+RSAVE(L-2,J,1))*DR1
40    CONTINUE
      D=1./(XJ*RL-XL*RJ)
      RR=XJ*D
      RX=-RJ*D
      XX=RL*D
      XR=-XL*D
      GO TO 46
45    CONTINUE
      X=XSAVE(L,J,1)
      R=RSAVE(L,J,1)
      CALL TRNX(X,R,XX,XR,CAPX,ZL)
      CALL TRNR(X,R,RR,RX,CAPR,RHUB)
46    CONTINUE
      XSAVE(L,J,2)=XR
      XSAVE(L,J,3)=XX
      RSAVE(L,J,2)=RR
      RSAVE(L,J,3)=RX
50    CONTINUE
100   CONTINUE
      RETURN
      END
      SUBROUTINE YYM(J)
      COMMON /XRTSAV/XSAVE(50,100,3),RSAVE(50,100,3)
      COMMON/TTSAVE/TTSAV(50,50,3),TRSAV(50,50),TXSAV(50,50),
     *   TPSAV(50,100),TSSAV(50,100)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/MESHSP/SPHALF,SKR,SKT
      REAL TJ(50),TL(50),TK(50)
      DT1=0.5/DTH
      NTH1=NTH-1
      DX1=0.5/DX
      DR1=0.5/DR
      JP=J+1
      JM=J-1
      DO 50 L=1,NR
      LP=L+1
      LM=L-1
      IF(J.EQ.1)GO TO 10
      IF(J.EQ.NX)GO TO 15
      XJ=(XSAVE(L,JP,1)-XSAVE(L,JM,1))*DX1
      RJ=(RSAVE(L,JP,1)-RSAVE(L,JM,1))*DX1
      DO 5 K=1,NTH
5     TJ(K)=(TTSAV(K,L, 3)-TTSAV(K,L, 1))*DX1
      GO TO 20
10    XJ=(-3.*XSAVE(L,J,1)+4.*XSAVE(L,J+1,1)-XSAVE(L,J+2,1))*DX1
      RJ=(-3.*RSAVE(L,J,1)+4.*RSAVE(L,J+1,1)-RSAVE(L,J+2,1))*DX1
      DO 13 K=1,NTH
13    TJ(K)=(-3.*TTSAV(K,L,1)+4.*TTSAV(K,L,  2)-TTSAV(K,L,  3))*DX1
      GO TO 20
15    XJ=( 3.*XSAVE(L,J,1)-4.*XSAVE(L,J-1,1)+XSAVE(L,J-2,1))*DX1
      RJ=( 3.*RSAVE(L,J,1)-4.*RSAVE(L,J-1,1)+RSAVE(L,J-2,1))*DX1
      DO 18 K=1,NTH
18    TJ(K)=(3.*TTSAV(K,L,3)-4.*TTSAV(K,L,2  )+TTSAV(K,L,  1))*DX1
20    IF(L.EQ.1) GO TO 30
      IF(L.EQ.NR)GO TO 35
      XL=(XSAVE(LP,J,1)-XSAVE(LM,J,1))*DR1
      RL=(RSAVE(LP,J,1)-RSAVE(LM,J,1))*DR1
      DO 25 K=1,NTH
25    TL(K)=(TTSAV(K,LP,2)-TTSAV(K,LM,2))*DR1
      GO TO 40
30    XL=(-3.*XSAVE(L,J,1)+4.*XSAVE(L+1,J,1)-XSAVE(L+2,J,1))*DR1
      RL=(-3.*RSAVE(L,J,1)+4.*RSAVE(L+1,J,1)-RSAVE(L+2,J,1))*DR1
      DO 33 K=1,NTH
33    TL(K)=(-3.*TTSAV(K,L,2)+4.*TTSAV(K,L+1,2)-TTSAV(K,L+2,2))*DR1
      GO TO 40
35    XL=(3.*XSAVE(L,J,1)-4.*XSAVE(L-1,J,1)+XSAVE(L-2,J,1))*DR1
      RL=(3.*RSAVE(L,J,1)-4.*RSAVE(L-1,J,1)+RSAVE(L-2,J,1))*DR1
      DO 38 K=1,NTH
38    TL(K)=(3.*TTSAV(K,L,2)-4.*TTSAV(K,L-1,2)+TTSAV(K,L-2,2))*DR1
40    CONTINUE
      DO 42 K=2,NTH1
42    TK(K)=(TTSAV(K+1,L,2)-TTSAV(K-1,L,2))*DT1
      K=1
      TK(1)=(-3.*TTSAV(K,L,2)+4.*TTSAV(K+1,L,2)-TTSAV(K+2,L,2))*DT1
      K=NTH
      TK(NTH)=(3*TTSAV(K,L,2)-4*TTSAV(K-1,L,2)+TTSAV(K-2,L,2))*DT1
      D=1./(XJ*RL-XL*RJ)
      RR=XJ*D
      RX=-RJ*D
      XX=RL*D
      XR=-XL*D
      XSAVE(L,J,2)=XR
      XSAVE(L,J,3)=XX
      RSAVE(L,J,2)=RR
      RSAVE(L,J,3)=RX
      DO 43 K=1,NTH
      D1=D/TK(K)
      TPSAV(K,L)=1./TK(K)
      TXSAV(K,L)=(TL(K)*RJ-TJ(K)*RL)*D1
43    TRSAV(K,L)=(XL*TJ(K)-TL(K)*XJ)*D1
50    CONTINUE
      RETURN
      END
      SUBROUTINE CALDER(NBLADE,IDMP)
       REAL X(50),R(50),XX(50),XR(50),RR(50),RX(50)
       REAL DTPDT(50),DTSDT(50)
       REAL TS(50),TP(50),DTPDS(50),DTSDS(50)
       REAL DTPDTP(50),DTSDTP(50),KS(2,50),KTAU(2,50)
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DZ,TAUX,TAUT,TAUR
      COMMON /BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
       REAL DTPDZ(50),DTSDZ(50),DTPDR(50),DTSDR(50)
       REAL CR(2,50),CT(2,50),CZ(2,50)
       REAL D2PDT(50),D2SDT(50),D2SDSP(50),D2PDSP(50),D2SDS(50)
       REAL D2PDS(50), DT2(2),DS2(2)
       DIMENSION DT(2),DS(2)
      REAL DTHZ(2),DTHR(2)
       REAL MOD, ALPT(20),ALPS(20)
       NR1=NR-1
      NEXTRA =0
      IF(IDMP.EQ.5)NEXTRA=1
      IF(IDMP.LT.0  .OR.  IDMP.EQ.6)NEXTRA=2
       TRAN=180.0/3.1416
       DO 100 J=N1,N2
      READ(2)X,R,XX,XR,RR,RX
       READ(4)TS,TP,DTPDT,DTSDT,D2PDT,D2SDT
       JP=J-N1+1
       DO 10 L=1,NR
      READ(8)DTPDTP,DTSDTP,D2PDSP,D2SDSP
       DTPDS(L)=DTPDTP(JP)
       DTSDS(L)=DTSDTP(JP)
       D2SDS(L)=D2SDSP(JP)
       D2PDS(L)=D2PDSP(JP)
10       CONTINUE
      REWIND 8
      DO 30 L=1,NR
C CALCULATE S AND TAU LINE ANGLES
       ALPS(L)=-RX(L)/RR(L)
       ALPS(L)=ATAN(ALPS(L))
       ALPT(L)=-XR(L)/XX(L)
       ALPT(L)=ATAN(ALPT(L))
C CALCULATE THE R & Z THETA DERIVATIVES
       COT=COS(ALPT(L))
       SIT=SIN(ALPT(L))
       COOS=COS(ALPS(L))
       SIS=SIN(ALPS(L))
       CSUM=COS(ALPT(L)+ALPS(L))
       DTPDZ(L)=(DTPDS(L)*COT-DTPDT(L)*SIS)/CSUM
       DTSDZ(L)=(DTSDS(L)*COT-DTSDT(L)*SIS)/CSUM
       DTPDR(L)=(-DTPDS(L)*SIT+DTPDT(L)*COOS)/CSUM
       DTSDR(L)=(-DTSDS(L)*SIT+DTSDT(L)*COOS)/CSUM
C CALCULATE THE R*THETA DERIVATIVES
      DELTH=TS(L)-TP(L)
      DTHZ(1)=DTPDZ(L)/DELTH
      DTHZ(2)=DTSDZ(L)/DELTH
      DTHR(1)=DTPDR(L)/DELTH
      DTHR(2)=DTSDR(L)/DELTH
       DT(1)=R(L)*DTPDT(L)
       DT(2)=R(L)*DTSDT(L)
       DS(1)=R(L)*DTPDS(L)
       DS(2)=R(L)*DTSDS(L)
C
       DT2(1)=R(L)*D2PDT(L)
       DT2(2)=R(L)*D2SDT(L)
       DS2(1)=R(L)*D2PDS(L)
       DS2(2)=R(L)*D2SDS(L)
C
       DO 40 I=1,2
C CALCULATE THE TAU VECTOR
       ROOT=SQRT(1+DT(I)**2)
       RT=COS(ALPT(L))/ROOT
       ZT=SIN(ALPT(L))/ROOT
       TT=DT(I)/ROOT
C CALCULATE THE S VECTOR
       ROOT=SQRT(1+DS(I)**2)
       RS=SIN(ALPS(L))/ROOT
       ZS=COS(ALPS(L))/ROOT
       TAS=DS(I)/ROOT
C CALCULATE THE NORMAL VECTOR
       RN=TT*ZS-ZT*TAS
       TN=ZT*RS-RT*ZS
       ZN=RT*TAS-TT*RS
C CALCULATE THE MODULUS
       MOD=SQRT(RN*RN+TN*TN+ZN*ZN)*(-1)**I
       IF(MOD.EQ.0.0)GO TO 40
C CALCULATE THE COSINE DIRECTORS
       CR(I,L)=RN/MOD
       CT(I,L)=TN/MOD
       CZ(I,L)=ZN/MOD
C CALCULATE THE CURVATURE
       KTAU(I,L)=DT2(I)/((1.0+DT(I)**2)**1.5)*(-1)**I
       KS(I,L)=DS2(I)/((1.0+DS(I)**2)**1.5)*(-1)**I
 40       CONTINUE
       IF(NEXTRA.NE.1)GO TO 30
       ALPS(L)=ALPS(L)*TRAN
       ALPT(L)=ALPT(L)*TRAN
       DTPDT(L)=TRAN*ATAN(DT(1))
       DTSDT(L)=TRAN*ATAN(DT(2))
       DTPDS(L)=TRAN*ATAN(DS(1))
       DTSDS(L)=TRAN*ATAN(DS(2))
 30       CONTINUE
      IF(NEXTRA.NE.1 .AND. NEXTRA.NE.2) GO TO 58
       IF(NEXTRA.NE.1)GO TO 51
      WRITE(6,600)J
600   FORMAT( 10X,' J= ',I5,' ALPT,ALPS IN DEG.')
       WRITE(6,201)(ALPT(L),L=2,NR1)
       WRITE(6,201)(ALPS(L),L=2,NR1)
      WRITE(6,601)
601   FORMAT(' PRESSURE SURFACE DTPDT DTPDS EXPRESSED IN DEG FROM',     -
     *   ' TAU OR S AXIS')
       WRITE(6,201)(DTPDT(L),L=2,NR1)
       WRITE(6,201)(DTPDS(L),L=2,NR1)
51    WRITE(6,602)
602   FORMAT(' CURVATURE IN TAU AND S DIRECTIONS ')
        WRITE(6,201)(KTAU(1,L),L=2,NR1)
        WRITE(6,201)(KS(1,L),L=2,NR1)
       IF(NEXTRA.NE.1)GO TO 52
      WRITE(6,603)
603   FORMAT(' DTSDT,DTSDS')
       WRITE(6,201)(DTSDT(L),L=2,NR1)
       WRITE(6,201)(DTSDS(L),L=2,NR1)
52    WRITE(6,602)
        WRITE(6,201)(KTAU(2,L),L=2,NR1)
        WRITE(6,201)(KS(2,L),L=2,NR1)
      IF(NEXTRA.NE.1) GO TO 58
      WRITE(6,604)
604   FORMAT(' DIRECTION COSINES')
       DO 50 I=1,2
       WRITE(6,201)(CR(I,L),L=2,NR1)
       WRITE(6,201)(CT(I,L),L=2,NR1)
       WRITE(6,201)(CZ(I,L),L=2,NR1)
 50       CONTINUE
 201       FORMAT(16F8.3)
58    CONTINUE
      WRITE(9)X,R,TS,TP,DTPDZ,DTSDZ,DTPDR,DTSDR
      WRITE(9)CR,CT,CZ
      WRITE(9)KS,KTAU
100       CONTINUE
      REWIND 9
      REWIND 2
      REWIND 4
      REWIND 8
      REWIND 3
       RETURN
       END
